On the critical behavior of the specific heat in glass-formers 
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We show numeric evidence that, at low enough temperatures, the potential energy density of a 
glass-forming liquid fluctuates over length scales much larger than the interaction range. We focus 
on the behavior of translationally invariant quantities. The growing correlation length is unveiled 
by studying the Finite Size effects. In the thermodynamic limit, the specific heat and the relaxation 
time diverge as a power law. Both features point towards the existence of a critical point in the 
metastable supercooled liquid phase. 
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The glass transition differs from standard phase transi- 
tions in that the equilibration time of glass- formers (poly- 
mers, supercooled liquids, colloids, granulars, super- 
conductors, ...) diverges without dramatic changes in 
their structural properties. Indeed elastic neutron scat- 
tering experiments show that the Static Structure Factor 
of polymers and supercooled liquids undergoes marginal 
modifications when approaching the glass transition. 
Yet, their equilibration time r, as obtained for instance 
from the a peak in the a.c. dielectric susceptibility, grows 
by many order of magnitude close to the transition[l|. 
Reconciliating the two faces of the glass transition is a 
major challenge for condensed matter physics. 

A general mechanism producing a divergence of the 
equilibration time of an homogeneous system at finite 
temperature is the divergence of a correlation length 
(critical slowing down^). Slowness is due to the fact 
that configurational changes need to propagate over in- 
creasingly large regions (the critical origin of the Mode 
Coupling singularity has been recently recognizedQ). 
Within this framework, a key problem is identifying the 
quantity with the largest spatial fluctuations close to the 
glass transition. Since equal-time correlation functions 
do not reveal statistical fluctuations over large length 
scales, it has been suggested that two-times correlation 
functions must be studied^, H, 0, 0, @- The goal is to 
extend beyond mean-field level the Mode Coupling view 
of the glass transition as a purely dynamic phenomenon. 

Strong space fluctuations of the relaxation prop- 
erties (dynamic heterogeneities) have been found 
experimentally 0, ^3 and numerically 0|- However, the 
size of the correlated domains is not larger than a few 
nanometers (a few angstroms, in simulations). 

Here we report numeric simulations for a fragile glass- 
forming liauidjl^. fl3| showing large scale fluctuations 
in the specific heat, an equal-time correlation function. 
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Thus, we claim that the dynamical features of the glass- 
transition are due to critical slowing down arising from 
a continuous phase transition suffered by the metastable 
liquid. The difficulty in recognizing it is due to the fact 
that standard scattering experiments are not devised to 
detect spatial fluctuations in the energy density. Yet, a 
large correlation-length can be studied through Finite- 
Size effect sll4l . Note that experiments in films |ll 5| and 
nanopores[l6J show that the glass transition changes in 
samples with one or more dimensions of nanometric scale, 
although it is difficult to disentangle Finite-Size Scaling 
from the effects of the interaction with the substrate. 
However, the specific-heat of toluene confined on pores of 
diameter 8.7 nm close to its glass temperature, is sig- 
nificantly smaller than for bulk toluene, which could sig- 
nal a correlation length well above the nanometric scale. 

We study a 50% mixture of particles interacting 
through the pair potential V a p(r) = e[(a a + ap)/r] 12 + 
C a p, where a, /3 = A, B, with a cutoff at r c = V3o~q. The 
choice (7b = 1-2(7,4 hampers crystallization. We impose 
(2cta) 3 + 2((7a + (7b) 3 + (2(7b) 3 = 4ctq where cto is the unit 
length. Constants C a p are chosen to ensure continuity 
at r c . The simulations are at constant volume, with par- 
ticle density fixed to (7q~ 3 and temperatures in the range 
[0.897T mc ,10.792T mc ], where T mc is the Mode Coupling 
temperature|3]. We use periodic boundary conditions on 
a box of size L (which discretizes momenta in units of 
q min = 2n/L) in systems with 512, 1024,2048 and 4096 
particles. From the Molecular Dynamics van Hove self- 
correlation function. 17], one has for argon parameters, 
(7 = 3.4A, e/k B = 120K and T mc = 26.4K. We shall 
obtain equilibrium data below T mc . 

We modify the Grigera-Parisi swap algorithm^] to 
make it local, in order to keep the algorithm in the 
dynamic Universality Class Q of standard Monte Carlo 
(MC). The elementary MC step is either (with proba- 
bility p) a single-particle displacement attempt or (with 
probability 1—p) an attempt to swap particles. The swap 
is made by picking a particle at random and trying to in- 
terchange its position with that of a particle of opposite 
type, chosen at random among those at distance smaller 
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than 0.6r c . Detailed balance requires that the Metropolis 
test include not only the energy variation but the change 
in the number of neighbors. The swap acceptance is inde- 
pendent of system size and grows from 0.74% at 0.9T mc 
up to 6% at 2T mc . In this work we use p = 0.5 (named 
local swap from here on) and p = 1.0 (named standard 
MC). The time unit to is N/p elementary steps. 

Although time correlators differ for different dynamics, 
when studying static quantities (e.g. specific heat) the 
choice of Monte Carlo dynamics is a matter of practical 
convenience. We study local swap time correlators only 
in order to explore static properties. Yet, the asymptotic 
decay of time correlators of two dynamics belonging to 
a single dynamic Universality Class is given by the same 
function (up to a rescaling) of the correlation length^. 
Since standard and swap MC are both local and share 
the same conservation laws, we expect the ratio of their 
autocorrelation times, Eq.Q, to be essentially constant, 
even if their time correlators differ at intermediate times. 

We assume that equilibrium behavior can be meaning- 
fully studied in a metastable phase^ij](i.e. the equilibra- 
tion time for the metastable liquid phase is much smaller 
than the crystallization time). At the lowest tempera- 
tures we simulate up to 400 independent MC runs. In 
the analysis we only consider histories longer than 100 
exponential autocorrelation times (see below). We use 
the Jack-Knife method to estimate errors |l4j. 

The study of a stochastic classical dynamics is math- 
ematically equivalent to the diagonalization of a quan- 
tum Hamiltonian[l9| . related to the Fokker-Planck (FP) 
operator, whose eigenvalues yield the relaxation time 
scales. Translational invariance of the FP operator in- 
duces a classification of physical quantities according to 
their wave vector, Q. The small Q density fluctuations 
are ruled by the hydrodynamic law stemming from mass 
conservation t(Q) oc Q _2 [2(J, which become singular at 
Q = 0. Yet, the dynamics of the translationally invari- 
ant quantities (Q = sector) may be studied for non 
conserved quantities such as the potential energy. 

A generic observable O is said to belong to the Q- 
sector if it transform as O — > e 1( ^' A under an uniform 
displacement of the coordinates, t; l — + A. Well known 
example are the Fourier transforms (Vjy = V(rk — 
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Note that p(0) is the (conserved) particle density, while 
p e (0) is the (non-conserved) potential energy density, 
which we call e (the internal energy is then |fc B T+(e)). 

Multiplying densities with wave- vectors q and —q 
yields translationally invariant observables: 

S{q) = p{q)p(-q) , S e {q) = p e {q)p e {-q) . (2) 

The mean value of S(q) is the static structure factor S(q), 
while (<S e (0)) is related to the constant-volume specific 




FIG. 1: Correlators for two translationally invariant quanti- 
ties: the potential energy density, e, and the square of the den- 
sity fluctuation, <S(g m in), at the minimal momentum allowed 
by the boundary conditions (1024 particles, below T mc ). The 
local swap algorithm decorrelates faster than standard MC. 



heat, Cy as T 2 C V = N{a ) 6 [(<S e (0)) - (e) 2 ] . For every 
observable, O, we consider the normalized time correlator 



C (t) = «O(0)O(i)> - (0) 2 )/((0 2 ) (O) 2 ) . (3) 

The Theory of Critical Phenomena^ suggest that, at 
very long times, time correlators decay exponentially |26(. 
We consider two autocorrelation timesji]]], the integrated 
time T* m t,0 an d the exponential time T cxp o ■ 



Tint.O 



C (t)dt, C {t) 



-t/r exp ,o 



(4) 



with different meanings. The Tj nti o of Q = quantities 
is related to transport properties (e.g. viscosityds T lrA for 
some components of the energy-stress tensor 20] ). The 
T C xp,o is the longest characteristic time and depends only 
on the Q-sector to which O belongs. 

In Fig. ^ we show the time correlator for the poten- 
tial energy density and <S(<f m in)- While at short times 
standard MC and local swap dynamics yield the same 
correlators, the swap does not present the Mode Cou- 
pling plateau (hence significance of T mc for swap MC is 
unclear), allowing a better study of the long time de- 
cay. Following this decay is difficult. The numerical (or 
experimental) effort to obtain the correlator C with pre- 
scribed accuracy grows as C~ 2 when C — > 0. In fact, 
most of previous work was confined to C > 0.1, while we 
are able to explore the range 0.01 < C < 0.1 (Fig.|2J. 

In Fig. [2] we show a semilogarithmic plot of the time 
correlators versus time. For long times, a straight line of 
slope — l/r cxp should be found. In Fig.^a) we show that 
T eX p obtained from e grows fast near the Mode Coupling 
temperature. The lower panels of Fig. [21 confirm that 
T cxp is actually a property of the Q = sector. When 
our statistics are good enough to follow the correlator for 
three decades, the estimate for the exponential time does 
not depend on the chosen correlator. At low temperature 
one must be content with choosing the most convenient 
observable to extract r cxp . Every observable within a 
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FIG. 2: Correlators for TV = 1024 particles (TV = 2048 at 
T/Tmc ~ 0.897) with the swap algorithm. Panel (a): tem- 
perature variation of the correlator of e, close to T mc . Panels 
(b)-(d): correlators of several Q = quantities. Asymptoti- 
cally, three parallel straight lines should be found, with slope 
— 1/Yexp- At the highest temperature, (b), this is clearly ob- 
served, as well as in (c). In (d) we do not have enough in- 
dependent measurements to see clearly this common slope. 
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i. reduced temperature (line is a fit 
of the swap data to a critical divergence). Note Universality 
with respect to changes in the dynamics (the two r cxp differ 
in a constant factor). With standard MC one is restricted to 
T > 0.973T mc . (b): ratio r int /r cxp for swap dynamics. The 
participation of <S(g m m) in the slowest mode decreases upon 
approaching T c , while the share of e slightly increases. We 
use TV = 1024 but for T/T mc = 0.897, where TV = 2048. 



symmetry class has some overlap with the slowest mode 
of that class. The one which shows the smaller slope at 
shorter times is closer to the slowest mode, allowing more 
precise studies of r exp . 

The growth of T exp is shown in Fig. OJa). Close to a 
critical point, one expects a power law divergence of T exp . 
This hypothesis accounts for our data, with an estimate 
for the critical temperature 0.84T mc < T c < 0.89T mc . 

Universality requires that T oxp oc A(T~T C )~ ZU , with T c 
and the combination of critical exponents zv independent 
of the dynamics. Interestingly enough, in Fig. E^a) we 
find that the ratio of r exp for the swap and the standard 
MC dynamics is constant within errors, in spite of the 
enormous difference on their short-time behaviors. 

A crucial issue is identifying the quantity with largest 
spatial fluctuations. For local dynamics, this quantity 
is closely related to the slowest mode (whose relaxation 
is purely exponential). Closeness can be quantified by 
the ratio of the integrated autocorrelation time to the 
exponential one, Eq.QJ. The challenge is to find the ob- 
servable that maximizes this ratio in the Q = sector 
(Tmt/7cxp = 1 only for the slowest mode). In FigOJb) 
we show the ratio Ti n t/V exp for the potential energy and 
for iS(g m in) , obtained with the swap dynamics. We ob- 
serve that S(q m m) seems a very good candidate at high 
temperatures but the ratio T in t / r exp sinks near T mc . For- 
tunately, the potential energy displays a modest but con- 
stant ratio (even slightly increasing at low temperatures). 
This is a strong indication that critical behavior can be 
investigated in the fluctuations of the potential energy, 
namely the specific heat, and that a diverging correla- 



tion length would show up in four-particle correlators. 

For the local dynamics studied here, diverging auto- 
correlation times suggest diverging correlation lengths, 
that should show up in static properties as well. Since 
the order parameter is unknown, a direct measure of 
the correlation length is difficult, but one may detect 
it indirectly through Finite-Size Scaling^]. The time 
correlators have taught us that energy fluctuations are 
promising candidates. We studied, Fig. 0Ja), the spe- 
cific heat dependence on the size of the simulation box, 
L. Unfortunately, small systems crystallize quickly be- 
low 2.13T mc . On the other hand, the metastable liquid 
can be studied with TV = 512 particles, using local swap, 
down to T = 0.897T mc . Up to T = 0.921T mc no finite- 
size effects are detected. However, for T = 0.897T mc , a 
noticeable growth of the specific-heat with L is found 
up to L — 12(To (~ 4nm for argon parameters). This 
length is comparable with the experimental domain size 
for cooperative dynamics^ 0], and well above previous 
simulations [llj. The r oxp show a similar effect, Fig.^Jb). 

The next step is the study of critical behavior in the 
infinite- volume specific heat itself, displayed in FigEJc). 
Generally speaking, critical divergences for the specific 
heat are difficult to study numerically due to the pres- 
ence of a large non-critical background (see e.g. Ref. 22]). 
Fortunately, in our case the background is given by the 
Rosenfeld-Tarazona law : 23], T 2 C V oc T 8 / 5 , which should 
be followed by a non-critical dense fluid at low temper- 
atures. We have checked that from T = 2T mc to be- 
yond 10T mc the T 8 / 5 law is extremely accurate. Interest- 
ingly, at lower temperatures (where the agreement with 
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FIG. 4: Panel (a): unlike internal energy (not shown), the 
specific heat is a growing function of the simulation box size, 
until the system becomes much larger than the correlation 
length^. This is visible for L < I2a at T = 0.897T mc , but 
not at T = 0.921T mc . The correlation length grows quickly 
close to T c , which is confirmed by the time correlator for the 
potential energy (b). Panel (c): infinite volume specific-heat 
vs. temperature. The dashed line is the y 8 / 5 Rosenfeld- 
Tarazona law. Deviations from the T 8//5 law can be fitted by 
a critical divergence (full line in (d)) yielding T c = 0.86T mc . 



the T 8 / 5 law should be still better in absence of criti- 
cality) deviations start to be significant. They are well 
described by a power law divergence oc (T — T c )~ a . To 
estimate errors in T c and a is difficult. Excluding the two 
extremal points in FigQ] (c), the fit yields a — 0.9 and 
T c = 0.86T mc . Including any of the extremal points in 
the fit, a slightly grows. Taking into account the bound 
a < 1 imposed by the continuity of the internal energy, 
we estimate 0.9 < a < 1.0 and 0.84T mc <T C < 0.86T mc . 

In summary, we studied the equilibrium static and dy- 
namic properties of a model of a fragile supercooled liq- 
uid, with emphasis on the dynamics of translationally in- 
variant quantities. We claim that critical slowing down is 
behind the structural arrest of glass formers, which have 
then universal properties. We find that time correlators 
with a complicated structure relax exponentially at very 
long times. The study of the time correlators is a pow- 
erful, unprejudiced method of identifying the physical 
quantities suffering fluctuations over the largest length 
scale. The potential energy, rather than density fluctu- 
ations, emerges as the candidate for the study of this 
critical phenomenon. Experiments measure the constant 
pressure specific-heat, that is obtained from the constant 
volume one by adding a term that is smooth in the ab- 
sence of critical density fluctuations (not found here). 
The critical temperature obtained from the divergence 
of the specific heat and of the autocorrelation times lies 
in the range T/T mc = 0.83—0.88 . The study of the dy- 
namics of translationally invariant quantities appears as 



a challenge to experimentalists. While measurements of 
the frequency dependence of the specific heat [3> 113 are 
an appealing possibility to estimate the potential energy 
relaxation time, the correlation-length could be studied 
by Finite-Size Scaling of the specific-heat and of relax- 
ation times in films 15] or in larg er pores than previously 
used to confine glass- formers jlq . 
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